#include "dEdxGausN.h"
#include "SpectraUtils.h"

TString GetYear(const TString fin)
{
  TString year;
  if(fin.Contains("CEN0"))
    year = "CEN0";
  else if(fin.Contains("CEN1"))
    year = "CEN1";
  else if(fin.Contains("CEN2"))
    year = "CEN2";
  else if(fin.Contains("CEN3"))
    year = "CEN3";
  else if(fin.Contains("CEN4"))
    year = "CEN4";
  else if(fin.Contains("CEN5"))
    year = "CEN5";
  else if(fin.Contains("10d") || fin.Contains("7000"))
    year = "10d";
  else if(fin.Contains("11a") || fin.Contains("2760")){
    year = "11a";
  }
  else{
    printf("bad year! %s\n", fin.Data()); exit(1);
  }

  return year;
}

void bwSpectra(TString fin, const Int_t ptmOpt, const Bool_t kfixprofile)
{

  Double_t ptmin[]={0,0,0};
  Double_t ptmax[]={100,100,100};
  if(ptmOpt==1){
    //official range
    ptmin[0]=0.3;
    ptmin[1]=0.5;
    ptmin[2]=0.2;

    ptmax[0] = 3;//proton
    ptmax[1] = 1;//pion
    ptmax[2] = 1.5; //kaon
  }
  else if(ptmOpt==2){
    //same range, different sample
    ptmin[0]=0.6;
    ptmin[1]=0.6;
    ptmin[2]=0.6;

    ptmax[0] = 3;//proton
    ptmax[1] = 1;//pion
    ptmax[2] = 1.5; //kaon
  }
  else if(ptmOpt==3){
    ptmax[0] = 3;//proton
    ptmax[1] = 1;//pion
    ptmax[2] = 1.5; //kaon
    /*
    ptmax[0] = 4;//proton
    ptmax[1] = 1.5;//pion
    ptmax[2] = 3; //kaon
    */
  }


  for(Int_t itype = 0; itype<3; itype++){
    printf("please check itype %d ptmin %e ptmax %e\n", itype, ptmin[itype], ptmax[itype]);
  }

  TFile * filin = new TFile(fin);

  TGraphAsymmErrors * grsys[]={0x0,0x0,0x0};
  TGraphAsymmErrors * grstat[]={0x0,0x0,0x0};
  Double_t mass[]={-999,-999,-999};
  TObjArray * arr = new TObjArray;

  const Int_t ntype = 3;
  for(Int_t itype = 0; itype<ntype; itype++){
    grsys[itype] = (TGraphAsymmErrors*)filin->Get(Form("grsys%d",itype));
    grstat[itype] = (TGraphAsymmErrors*)filin->Get(Form("grstat%d",itype));
    if(!grsys[itype] || !grstat[itype]){
      printf("no gr sys or stat %d\n", itype);filin->ls(); exit(1);
    }

    mass[itype]=dEdxGausN::Mass(itype);
    printf("mass %d %e\n", itype, mass[itype]);
  }

  for(Int_t itype = 0; itype<ntype; itype++){
    TGraphAsymmErrors * grcomb = new TGraphAsymmErrors;
    grcomb->SetName(Form("grcomb%d",itype));

    for(Int_t ip=0; ip<grsys[itype]->GetN(); ip++){
      const Double_t xx = grstat[itype]->GetX()[ip];
      const Double_t xx2 = grsys[itype]->GetX()[ip];
      if(fabs(xx-xx2)>EPSILON){
        printf("wrong xx %d %d: %e %e\n", itype, ip, xx, xx2); exit(1);
      }

      //do the range selection only in _FCN
      /*
      if(xx<ptmin[itype] || xx>ptmax[itype]){
        continue;
      }
      */

      const Double_t exlow = grstat[itype]->GetEXlow()[ip];
      const Double_t exlow2 = grsys[itype]->GetEXlow()[ip];
      if(fabs(exlow-exlow2)>EPSILON){
        printf("wrong exlow %d %d: %e %e\n", itype, ip, exlow, exlow2); exit(1);
      }

      const Double_t exhigh = grstat[itype]->GetEXhigh()[ip];
      const Double_t exhigh2 = grsys[itype]->GetEXhigh()[ip];
      if(fabs(exhigh-exhigh2)>EPSILON){
        printf("wrong exhigh %d %d: %e %e\n", itype, ip, exhigh, exhigh2); exit(1);
      }

      //==================
      //restore to dN instead of
      //dN/2/pi/pt which is read in
      const Double_t ff = 2*TMath::Pi()*xx;

      const Double_t yy = grstat[itype]->GetY()[ip]*ff;
      const Double_t yy2 = grsys[itype]->GetY()[ip]*ff;
      if(fabs(yy-yy2)>EPSILON){
        printf("wrong yy %d %d: %e %e\n", itype, ip, yy, yy2); exit(1);
      }

      const Double_t eylow = grstat[itype]->GetEYlow()[ip]*ff;
      const Double_t eylow2 = grsys[itype]->GetEYlow()[ip]*ff;
      const Double_t eyhigh = grstat[itype]->GetEYhigh()[ip]*ff;
      const Double_t eyhigh2 = grsys[itype]->GetEYhigh()[ip]*ff;
   
      //PWG without sys or stat then it is their quadratic sum!
      const Double_t eylowC = TMath::Sqrt(TMath::Power(eylow,2)+TMath::Power(eylow2,2));
      const Double_t eyhighC = TMath::Sqrt(TMath::Power(eyhigh,2)+TMath::Power(eyhigh2,2));

      const Int_t tmpnp = grcomb->GetN();
      grcomb->SetPoint(tmpnp, xx, yy);
      grcomb->SetPointError(tmpnp, exlow, exhigh, eylowC, eyhighC);
      printf("setting grcomb %d %d: %e %e , %e %e %e %e\n", itype, tmpnp, xx, yy, exlow, exhigh, eylowC, eyhighC);
    }

    arr->Add(grcomb);
  }

  TList *lout = new TList;
  lout->SetName("lout");

  const TString year=GetYear(fin);

  TList * ll = BGBlastWave_GlobalFit(year, arr, mass, ptmin, ptmax, lout, kfixprofile?1:0.7, kfixprofile);
  ll->ls();
  lout->ls();

  while(fin.Contains("/")){
    fin=fin(fin.First('/')+1,fin.Length());
  }

  const TString pwd=gSystem->pwd();
  TString tag="bw";
  if(pwd.Contains("POS"))
    tag+="_POS";
  else if(pwd.Contains("NEG"))
    tag+="_NEG";

  tag+=Form("_ptmOpt%d_fixPr%d_%s_PWG%d.root", ptmOpt, kfixprofile, year.Data(), pwd.Contains("PWG"));


  TFile * filout = new TFile(tag,"recreate");

  tag.ReplaceAll(".root","");
  ll->Write(tag, TObject::kSingleKey);
  lout->Write();
  filout->Save();
  filout->Close();
}


int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=4){
    printf("argc!=2\n"); return 1;
  }

  bwSpectra(argv[1], atoi(argv[2]), atoi(argv[3]));

  return 0;
}
